#Helper function for doing bootstraps for implementation results. 
cluster_bootstrap<-function(base,n,amn,bin){
  var_names<-unlist(strsplit(
    unlist(strsplit(colnames(base$model),
                    split="as.factor(", fixed=T)), split=")", fixed=T))
  model_data<-na.omit(base$data[!rownames(data) %in% na.action(base),
                        c(var_names, "pam_caseid")])
  pop<-unique(base$data$pam_caseid)
  coef<-matrix(NA, ncol=3)
  i<-0
  model<-NA
  while(length(na.omit(coef[,2]))<n){
    i<-i+1
    #Make within loop holder for effects
    effects<-matrix(NA, ncol=3)
    sample<-sample(pop, replace = T, size = length(pop))
    try(map<-sapply(sample, function(x) which(model_data[,"pam_caseid"]==x)))
    try(boot.sample<-data.frame(model_data[unlist(map),]))
    #Check for data that's not full rank
    if(any(apply(model_data,2,sd)==0)) next
    #Full Controls
    try(model<-brglm(formula(base), data=boot.sample))
    if(any(is.na(model))) next
    if(any(is.na(model$coefficients))) next
    if(model$converged==F) next
    atlist<-list(dummy=c(0,1))
    names(atlist)<-bin
    effects[,1:2]<-summary(margins(model, variable=paste(amn), at=atlist))[,3]
    effects[,3]<-effects[,2]-effects[,1]
    #cat("\r", i, "of", n) 
    coef<-rbind(coef, effects)
    model<-NA
    boot.sample<-NA
  }
return(coef)
}
